PL-TR-97-2031 


SURFACE  WAVE  WAVEFORM  INVERSION 
USING  GDSF  THEORY 


R.  B.  Herrmann 


Department  of  Earth  and  Atmospheric  Sciences 
St.  Louis  University 
3507  Laclede  Avenue 
St.  Louis,  MO  63103 


28  June  1996 


Scientific  Report  No.  1 


Approved  for  public  release;  distribution  unlimited 


DEPARTMENT  OF  ENERGY 
Office  of  Non-Proliferation 
and  National  Security 
WASHINGTON,  DC  20585 


19980819  091 


PHILLIPS  LABORATORY 
Directorate  of  Geophysics 
AIR  FORCE  MATERIEL  COMMAND 
HANSCOM  AFB,  MA  01731-3010 


SPONSORED  BY 
Department  of  Energy 

OfBce  of  Non-Proliferation  and  National  Security 

MONITORED  BY 
Phillips  Laboratory 
CONTRACT  No.  F19628-95-K-0019 

The  views  and  concliisions  contained  in  this  document  are  those  of  the  authors  and 
should  not  be  interpreted  as  representing  the  official  policies,  either  express  or  implied,  of 
the  Air  Force  or  U.S.  Government. 

This  technical  report  has  been  reviewed  and  is  approved  for  publication. 


JAMES  F.  LEWKOWIC^ 
Alternate  Contract  Mana^r 
Earth  Sciences  Division 


IS  F.  LEWKOWICZ 
Director 

Earth  Sciences  Division 


This  report  has  been  reviewed  by  the  ESD  Public  Affairs  Office  (PA)  and  is  releasable  to 
the  National  Technical  Information  Service  (NTIS). 


Qualified  requestors  may  obtain  copies  fi:om  the  Defense  Technical  Information  Center. 
All  others  should  apply  to  the  National  Technical  Information  Service. 


If  your  address  has  changed,  or  you  Avish  to  be  removed  from  the  mailing  list,  or  if  the 
addressee  is  no  longer  employed  by  your  organization,  please  notify  PL/IM,  29  Randolph 
Road,  Hanscom  AFB,  MA  01731-3010.  This  will  assist  us  in  maintaining  a  current 
mailing  list. 


Do  not  return  copies  of  this  report  unless  contractual  obligations  or  notices  on  a  specific 
document  requires  that  it  be  returned. 


Form  Approved 
0MB  No.  0704-0188 


REPORT  DOCUMENTATION  PAGE 

Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources, 
gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this 
collection  of  information,  including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson 
Davis  Highway.  Suite  1204,  Arlington,  VA  22202-4302,  and  to  the  Office  of  Management  and  Budget,  Paperwork  Reduction  Project  (0704-0188),  Washington,  DC  20503. 

AGENCY  USE  ONLY  (Leave  [  2.  REPORT  DATE  3.  REPORT  TYPE  AND  DATES  COVERED 

_  28  June  1996  Scientific  Report  No,  1 


4.  TITLE  AND  SUBTITLE 

Surface  Wave  Waveform  Inversion  Using  GDSF  Theory 

5.  FUNDING  NUMBERS 

Contract  No: 

F19628-95-K-0019 

PE  69120H 

6.  AUTHORS 

R.B.  Herrmann 

PR  DENN 

TA  GM 

WU  AJ 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

8.  PERFORMING  ORGANIZATION  REPORT 

Department  of  Earth  and  Atmospheric  Sciences 

NUMBER 

Saint  Louis  University 

3507  Laclede  Avenue 

St.  Louis,  MO  63103 

9.  SPONSORING  / MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

10.  SPONSORING /MONITORING  AGENCY 

Phillips  Laboratory 

REPORT  NUMBER 

29  Randolph  Road 

PL-TR-97-2031 

Hanscom  AFB,  MA  01731-3010 

Contract  Manager:  Delaine  Reiter/GPE 

11.  SUPPLEMENTARY  NOTES 


This  research  was  sponsored  by  the  Department  of  Energy,  Office  of 
Non-Proliferation  &  National  Security,  Washington,  DC  20585 

12a  DISTRIBUTION /AVAILABILITY  STATEMENT  12b.  DISTRIBUTION  CODE 

Approved  for  public  release;  distribution  unlimited 


13.  ABSTRACT  (Maximum  200  words) 

The  use  of  regionally  recorded  Rayleigh  waves  to  estimate  isotropic  moment  and 
psi-infinity  for  the  purpose  of  yield  estimation  is  investigated.  These  seismological 
parameters  and  their  variability  are  consistent  with  other  investigations.  Psi- 
infinity  shows  less  variability  than  isotropic  moment  as  a  yield  estimator.  Preliminary 
results  on  wave  propagation  in  the  Korean  peninsula  are  also  presented. 


14.  SUBJECT  TERMS  15.  NUMBER  OF  PAGES 

Isotropic  moment  -  yield,  nuclear  explosion,  chemical  _ 

explosion,  seismic  attenuation  16.  PRICE  CODE 

17.  SECURITY  CLASSIFICATION  I  18.  SECURITY  CLASSIFICATION  19.  SECURITY  CLASSIFICATION  1  20  LIMITATION  OF  ABSTRACT 
OF  REPORT  OF  THIS  PAGE  OF  ABSTRACT 

_ Unclassified _ Unclassified _ Unclassified 

NSN  ^4001-2806500 


Standard  Form  298  (Rev.  2-89) 
Prescribed  by  ANSI  Std.  Z39-1 
298-102 


SECTION 


TABLE  OF  CONTENTS 


PAGE 


Summary  ^ 

SURFACE  WAVE  WAVEFORM  INVERSION  USING  GSDF  THEORY  1 


1. 

Introduction 

1 

2. 

Theory  of  Generalized  Seismological  Data  Functionals 

1 

3. 

Structure  Inversion 

7 

4. 

Forming  Inversion  Kernel 

16 

5. 

Synthetic  Test 

20 

6. 

Real  Data  Test 

21 

7. 

Inversion  Procedure 

23 

8. 

Inversion  Results 

29 

9. 

Conclxision 

33 

10. 

References 

34 

iii 


This  page  left  blank 


IV 


Summary 

This  report  consists  of  a  review  of  the  interpretation  of  Generalized  Seis- 
mological  Data  Functionals  (GSDF)  as  a  means  of  refining  earth  structure 
through  waveform  modeling.  An  inversion  procedure  is  developed,  and  exam¬ 
ples  are  shown  for  synthetic  and  observed  data. 

The  focus  has  been  on  the  algorithm  implementation,  and  more  work 
must  be  done  to  imderstand  the  limitations  and  applicability  of  the  technique 
with  respect  to  real  data.  The  effects  poorly  constrained  source  depth  and 
mechanism  as  well  as  noisy  data  remain  to  be  investigated.  A  useful  tool  has 
been  implemented,  though. 
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SURFACE  WAVE  WAVEFORM  INVERSION  USING  GSDF  THEORY 

Tao-Ming  Chang  and  Robert  B.  Herrmann 

1.  Introduction 

In  this  report,  we  extend  the  Generalized  Seismological  Data  Function¬ 
als  (GSDF)  theory  to  the  inversion  of  broadband  waveforms  modeled  as  a 
superposition  of  sinface-wave  modes.  We  demonstrate  the  utility  of  this 
inversion  algorithm  by  conducting  a  simple  S3nithetic  test,  and  then  apply  the 
algorithm  to  real  data  to  see  what  new  knowledge  about  earth  structure  can 
be  obtained  from  this  new  inversion  algorithm. 

2.  Theory  of  Generalized  Seismological  Data  Functionals 

Gee  and  Jordan  (1992)  introduced  the  theory  of  Generalized  Seismologi¬ 
cal  Data  Functionals  (GSDF)  as  a  means  of  characterizing  differences 
between  observed  and  predicted  seismograms  in  terms  of  time  delays  with 
respect  to  the  earth  model  used.  The  importance  of  the  Gee  and  Jordan 
(1992)  paper  is  that  it  defines  the  framework  for  working  with  multimode  sig¬ 
nals. 

To  use  the  GSDF  approach,  we  need  to  construct  a  synthetic  seismogram 
is),  an  isolation  function  {f)  and  single-mode  seismograms.  An  isolation 
function  is  a  sum  of  single-mode  seismograms  which  may  represent  the  domi¬ 
nant  part  of  the  observed  seismogram  (s).  Using  a  cross-correlation  tech¬ 
nique,  we  can  quantify  the  similarity  observed  and  S3mthetic  seismograms. 
For  a  model  which  does  not  significantly  deviate  from  actual  earth  structmre, 
the  peak  of  cross-correlagrams  will  be  located  near  zero  lag-time  for  all  fre¬ 
quency  bands  and  for  all  windowed  segments  of  the  seismograms.  The  cross- 
correlagrams  reflect  the  degree  that  the  model  used  to  create  synthetic  seis¬ 
mogram,  isolation  function  and  single-mode  seismograms  is  different  from 
real  earth  structure.  To  utilize  this  information,  we  have  to  extract  informa¬ 
tion  from  different  frequency  ranges  and  interpret  it  in  term  of  differences  in 
earth  structure. 
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Gee  and  Jordan  (1992)  call  the  extracted  information  Generalized  Seis- 
mological  Data  Fimctionals  (GSDF).  In  the  following,  we  give  a  brief  review 
of  how  the  GSDF  is  defined,  how  the  GSDF  is  related  to  familiar  physical 
quantities,  and  how  we  extend  this  theory  to  invert  for  Earth  structure  Equa¬ 
tions  directly  adapted  fi”om  Gee  and  Jordan  (1992)  will  be  cited  as  (GJ.#), 
where  #  indicates  the  original  equation  number  in  their  paper. 

First,  we  review  the  definition  and  computation  of  a  GSDF.  We  can  use  a 
Gaussian  wavelet  to  approximate  the  filtered  and  windowed  cross- 
correlagrams  at  a  specified  frequency.  Specifically,  the  windowed,  filtered 
cross-correlagram  of  an  isolation  function  with  a  seismogram  (either  observed 
or  synthetic)  is  modeled  using  a  five-parameter  Gaussian  wavelet: 

FiWCfsit)  =  g(t)  =  AGa[cTs(f  -  tg)]  cos[<ys(i  -  tp)]  (GJ.5) 

YiVfCfsit)  ~  g{t)  =  AGa[<Ts(i  -  f^)]  cos[tys(^  -  tp)]  (GJ.8) 

F^WC^(^)  =  G8i[dft]  cos[& ft] 

where 

Fj  is  a  Gaussian  frequency  content  filter  with  center  frequency  coi  and 
half-bandwidth  cTj, 

^  ~  +  0- 01[o’^(^  ~  • -j  is  a  temporal  window  with 

half-bandwidth  c7„,and  centered  at  usually  t^  =  t^,  where  t^  is  the  lag¬ 
time  of  the  peak  of  FjWC^y, 

=  fit)  sit),  0  represents  cross-correlation,  is  the  cross-correlation 
with  the  observed  time  history, 

Cfsit)  =  fit)  0  sit),  is  the  cross-correlation  with  the  synthetic  time  his¬ 
tory, 

Ga[3c]  =  exp(-  —),  Ga  is  the  Gaussian  function, 

A  is  the  amplitude  of  Gaussian  envelope. 
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a  is  the  half-bandwidth  of  envelope  spectrum, 

CO  is  the  angular  frequency  of  oscillating  wavelet, 
tg  is  the  envelope  group  delay  from  zero  lag-time, 
tp  is  the  wavelet  phase  delay  from  zero  lag-time, 

the  subscript  s  denotes  observed  seismogram,  subscript  s  combined  with 
-  denotes  sjmthetic  seismogram,  and  the  subscript  f  with  ~  denotes  iso¬ 
lation  filter. 

From  (GJ.5)  and  (GJ.8),  Gee  and  Jordan  define  four  data  functionals  to  char¬ 
acterize  the  agreement  between  observed  and  predicted  time  histories; 


II 

1 

(GJ.9) 

~^g  ^ g 

(GJ.IO) 

1 

St„  -  -  -WnA-lnA] 

(GJ.ll) 

^ta=-^\.COs-^s\ 

(GJ.12) 

These  four  data  functionals  are  related  to  differential  phase  delay,  differen¬ 
tial  group  delay,  differences  in  logarithmic  amplitudes  and  the  differences  in 
center  frequencies.  These  four  data  functionals  are  defined  from  two  filtered, 
windowed  cross-correlagrams,  they  can  be  transformed  into  the  more  familiar 
quantities  such  as  phase  velocity,  group  velocity,  and  attenuation  corrections. 

We  will  briefly  describe  this  transformation.  We  note  that  all  filtered 
and  windowed  cross-correlagrams  will  be  normalized  by  scaling  FjWC^  to 
iinit  amplitude.  Also,  if  the  observed  and  synthetic  seismograms  are  com¬ 
posed  only  of  a  single  mode,  then  the  GSDF  are  easily  interpreted  in  terms  of 
differences  in  modeled  and  actual  phase  velocity,  group  velocity  and  Q.  For 
multi-mode  time  histories  their  interpretation  is  much  more  difficult. 

Before  relating  GSDF  to  physical  quantities,  we  must  know  the  roles 
these  quantities  play  in  wave  propagation.  For  a  known  instrument  response 
and  source,  the  difference  between  the  spectrum  of  isolation  filter 
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fia)  =  I{a)P{co)S{co)  and  its  corresponding  component  of  the  observed  seismo¬ 
grams  ficai)  =  I{g))P{(o)S{co)  is  called  "differential  propagation".  The  are 
related  by  the  equation 


=  D{(i))f{co)  (GJ.42) 

where  R  is  source-receiver  distance.  To  a  first  order  approximation,  the  dif¬ 
ferential  propagation  can  be  approximated  as 


D(C0)  =  exp  -COjdTqieOj)  -{co-  0)j)STa(cOj) 


(GJ.44) 


+exp  ^(OjSTpioj)  +  (eo-  o)  j)STq(o)j)^ 


The  propagation  effects  of  the  isolation  filter  fit)  and  the  corresponding 
waveform  f(t)  on  the  observed  seismogram  are  defined  in  four  equations 
which  are  described  in  terms  of  familiar  physical  quantities  such  as  total 
phase  delay,  total  group  delay,  and  attenuation  factor. 


SXpicOj)  =  XpicOj)  -  XpicOj) 

(GJ.45) 

STgiaOj)  =  TgicOj)  -  ig{(Oj) 

(GJ.46) 

Stqioj)  =  ^  fp{a>j)[Q-^  -  Q  -1] 

(GJ.47) 

Sxaicoj)  =  ^  fg((i}j)[Q-^  -  Q  -1] 

(GJ.48) 

where 

Tpicoj)  is  total  phase  delay  of  f(t)  at  an 
respect  to  the  original  time, 

fpicoj)  is  total  phase  delay  of  fit)  at  coj. 

arbitrary  frequency  oj  with 
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Tgiojj)  is  total  group  delay  of  fit)  at  (Oj, 
fgicoj)  is  total  group  delay  of  fit)  at  0j, 

Q  is  attenuation  factor  of  fit), 

Q  is  attenuation  factor  of  fit). 

This  means  that  whenever  we  are  able  to  measure  four  differential  quanti¬ 
ties  Stx  from  observed  and  synthetic  seismograms  at  a  certain  frequency  cdj, 
we  can  approximate  the  true  waveform  spectrum  behavior  fico)  at  any  arbi¬ 
trary  frequency  <y  to  a  first  order  precision. 

The  GSDF  St,,  based  on  the  filtered,  windowed  cross-correlation,  and  the  St^; 
is  related  to  the  time  domeiin  windowing  function.  Their  relationships  are 
given  by 


Stg  =  StgiSf)  +  (1  -  -  Stgi&f)] 

(GJ.56) 

Stp  *  Stpiwf)  +  (1  -  f,)i^l^)[t,  -  Stgi&f)] 

COf 

(GJ.57) 

Sta  ==  ^ISTai&f) 

(GJ.58) 

6tq~  S-tqi&f)  (1  ^f)(  )6T:ai&f) 

(Of 

(GJ.59) 

where 

is  the  frequency  from  Cff, 
is  a  window  width  factor, 

tc  is  the  lag-time  of  the  peak  of  cross-correlagram  FjWC^. 

There  is  an  additional  assumption  that  must  be  stated.  Due  to  the  difficulty 
in  relating  the  isolation  filter’s  corresponding  feature  in  observed  seismo¬ 
gram,  it  is  not  easy  to  evaluate  FjWC^.  However,  if  the  windowing  proce¬ 
dure  is  effective  in  isolating  fit)  from  other  phases  on  the  seismograms,  then 
F,WC^,=F,WC^. 

We  can  estimate  t^  from  an  isolated  waveform.  For  non-isolated  wave¬ 
forms,  the  Cff  can  not  be  simply  replaced  by  Cyj  and  this  is  a  more  general 
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situation.  For  non-isolated  waveforms,  the  Gaussian  wavelet  can  be  repre¬ 
sented  as  a  sum  of  interference  eifects  of  different  modes. 

¥iWCf,(t)  =  i  £  F,WC,„„(^)  (GJ.63) 

m-1 71=1 


FiWC^(i)=  f  X  FiWC^W 

m=l  n=l 

where 

=  [mode  m  off(t)]  gi  [mode  n  of 
=  [mode  m  off  if)]  iS  [  mode  n  of  fif)]. 

We  now  characterize  these  by  the  relation 

FiWC„,,,a)  »  A^,^GB.[df{t  -  ^y")]  cos[®,„„(^  -  tj")]  (GJ.68) 

where 


^mn=<^f  a  • 

In  this  case,  we  need  to  know  four  time  shifts  that  describe  the  deviations  of 
FjWC^s  from  FjWCy5r.  Assume  that  the  windowing  and  filtering  effectively 
suppress  the  bandwidth  variations,  so  that  ==  df.  The  four  time  shifts  are 
the  phase  delay  tp,  a  group  delay  tg,  and  two  amplitude  parameters 


t 


9  ~ 


-^InA 

(Df 


(GJ.60) 


<^f 

By  defining  the  following  notation, 

^mn  ~  2^^^  exp[— g  —t  g)]Ga[<Ty(^  g  —t  g)] 


(GJ.61) 


(GJ.73) 


-  ^f(tT  -  taXtp  -  tg)  (GJ.74) 

the  perturbation  expansions  can  be  simplified  by  using  a  set  iV  x  iV  matrices 

(C) 

mn  ~  Bmn  COS  ^mn  (GJ.75) 
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(S);„„  =  sin  (GJ  .76) 

-  t\)cos(Pmn  (GJ.77) 

(Sx)mn  =  ^mn®f  (C”  ”  t^)sin<l)mn  (GJ -78) 

Gee  and  Jordan  (1992)  gave  the  following  linearized  relationships  between 
the  GSDFs  and  computable  quantities  from  individual  mode  branches. 

=  1  •  C  •  +  1  •  S  ■  (GJ.84) 

=  - 1  •  S  •  ^tq  +  1  •  C  •  ^tp  (GJ.85) 

dta=-l-  (Ca  +  Sg)  •  ^tq  +  1  •  (Cg-Sa)  •  <Jtp  +  1  •  C  •  +  1  •  S  •  S%  (GJ.86) 

5tg=-l-  (Cg-Sa)  •  Jtq  -  1  -  (Ca  +  Sg)  •  Jtp  -  1  •  S  -  <5ta  +  1  •  c  •  (GJ.87) 


The  1  is  a  AT-dimension  vector,  each  element  of  which  equals  one.  The  is 
a  JV-vector  whose  n’th  component  is  dt'p.,  the  perturbation  corresponding  to 
the  quantity  x  of  n’th  mode.  In  the  following  section,  we  will  translate  this 
into  St^  which  directly  relate  to  seismogram  and  has  clear  physical  mean¬ 
ing. 

3.  Structure  Inversion 

To  apply  GSDF  theory  in  structure  inversion,  we  must  construct  the 
inversion  kernel,  G  for  our  inversion.  The  inverse  problem  can  be  simply 
expressed  in  the  form  : 

Jtx  =  Gx  • 

where 

X  indicates  one  of  {p,  g,  q,  a}. 

(Jtx  is  an  iV-vector  for  corresponding  to  N  modes, 

is  an  iV  X  ^  Frechet  kernel  matrix  for  structural  inverse  problem, 

(5m  is  a  model  correction  vector  for  k  unknowns. 

Keep  in  mind,  that  is  not  measurable  but  fortunately  GSDF  theory  pro¬ 
vides  a  way  to  compute  these  nonmeasurable  quantities  using  the  mode 
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interference  relationship.  Thus,  in  actual  application  of  GSDF  theory,  we 
must  relate  to  to  create  the  kernel  G^. 

From  (GJ.56-59),  we  can  transform  to  as 

Stg  =  SxgiSf)  +  (1  -  STgi^f)] 

Stp  =  Sxpidif)  +  a[tcl  -Sxgi&f)] 

St^  = 

^tq  =  5zq(^&f)- aSrJ^af) 

where  a  =  (1  -  Substituting  these  into  (GJ.84-87)  gives  the  fol- 

COf 

lowing  equations  which  relate  the  four  GSDFs  to  inversion  kernels  G^. 

Stp  =  1  •  C  •  la  +  1  •  ^C(Gp  -  aGg)  —  S(Gq  -  aGj,)j  <5m 

Stq  =  1  •  S  •  la  +  1  •  |^C(Gq  -  aGa)  +  S(Gp  -  aGg)J 
St^  =  l-  (Cg-Sa)  •  la  +  1  •  S  •  l(l-^f)i, 

+ 1  •  [-(C.  +  SjXG,  -  oG.)  +  (Cj  -  S.XGp  -  oGj)  +  If  CG.  +  |f  SG J 
=  - 1  ■  (Ca  +  Sg)  ■  la  +  1  •  C  •  1(1  -  ^l)t^ 

+  ij^— (Cg  —  Sa)(Gq  —  aGa)  ~  (Ca  +  Sg)(Gp  —  oGg)  —  SGa  +  ^fCGgj 

Now,  we  have  a  general  form  of  inversion  problem  for  earth  structure  in 
terms  of  the  GSDF’s  St^,  kernels  G^,  and  the  interference  effects  C,  S,  Cx,Sx. 

At  this  point,  it  is  possible  to  outline  this  structure  inversion  problem. 
At  a  particular  frequency,  according  to  GSDF  theory,  using  cross-correlation 
technique  we  can  obtain  foiir  measurements  as  our  "observations"  (GJ.9-12) 
and  evaluate  the  mode  interference,  which  we  will  use  with  G^  to  form  inver¬ 
sion  kernels.  Therefore,  we  perform  this  measuring  procedure  at  several  fre¬ 
quencies,  we  then  have  the  "data  vector"  and  "kernel  matrix"  ready  for 
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inversion.  In  the  following  sections,  we  will  show  how  we  create  the  kernels 
(GJ. 


Kernel  Gp 

Kernel  Gp  relates  phase  velocity  changes  to  model  perturbations. 

^Tp  =  Gp  ■ 

The  n’th  component  in  the  differential  phase  velocity  vector  5z^,  Szp^  is  define 
as 


R 


'nobs 


R 


c 


nsyn 


R 


R  R 


Cji  + 


££iL  +  (££5:)2 


where 


R 


R  dCn 
c„2  9m 


- 


R  is  soimce-receiver  distance, 

c„  is  phase  velocity  for  n’th  mode  of  synthetic  seismogram. 

In  the  matrix  form,  we  can  exactly  see  the  meaning  of  5Xp^  as  the  total  phase 
delay  for  mode  1  when  model  m  changes  into  m  +  ^m. 


R  dci 

R  dci 

R  dci 

- 

Sxp^ 

-  - 

c\  dmi 

R  dc2 

Cj  dm2 

R  dco 

c\  dm^ 

R  3co 

(#1  mode) 

'  Smi ' 

5  m2 

Nxl 

c|  dmi 

Cg  dm2 

c\  dm^ 

(#2  mode) 

Sm^ 

R  dcjf 
.  c%  dmi 

R  dc^ 
c%  dm2 

R  dc^ 
c%  dms 

■  •  •  (#iV  mode) 

Nxk 
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Kernel  Gg 


To  construct  kernel  Gg  which  relates  group  velocity  variation  to  model 
perturbations, 


STg  =  Gg  •  Sia. 


the  n’th  component  in  differential  group  velocity  vector  Stg,  dzg^  is  defined  as 


Stg  = 

U 


R 


R 


Jiobs 


u 


nsyn 


R 


R  R 


U,  +  SU,  £/„  U, 


U 


u. 


R 

U„ 


TT  '  V  T-T  / 


u. 


Ur. 


uj  ' 


R  dUr 


UJ  am 


^m 


where 

is  group  velocity  of  n’th  mode  synthetic  seismogram. 
In  matrix  form 


■  ■ 

^^g2 

R  aC7i 
Ul  ami 

R  dU2 
Ul  dmi 

R  BUi 
Ul  dm2 

R  BU2 
Ul  Bm2 

R  at/i 

Ul  Bms 

R  BU2 
Uldm, 

•  •  •  (#1  mode) 

•  •  •  (#2  mode) 

6mi ' 

5  m2 
Sm^ 

■  ^^gK  - 

Nxl 

R  BUn 
.  U%  9mi 

R  BUn 

u%  3^2 

R  Wt, 

V%  3^3 

■  ■  ■  {#N  mode) 

Nxk 

_  Sm^  _ 
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Kernel  G„ 


To  create  kernel  Gq  which  relates  attenuation  to  model  perturbations, 


<5rq  =  Gq-^m 


we  extend  (GJ.47)  into  a  general  matrix  form  and  reformulate  the  above 
equation  to  yield 


irka>) 


1-1 


dm 


(^m 


=  Gq  • 

Using  the  relationship  between  y  and  Q,  the  peirtial  derivative  of  y  with 
respect  to  model  perturbations  can  be  calculated. 

CO 


rn  = 


2C(i 


^Yn  ^  ^Q~n  _  Q-1  ^ 

dm  2co  dm  "  2co  ^  dm 
Therefore,  the  partial  derivative  of  can  be  rewritten  as 


dm 


^Yn  Q-1  9C0„ 

dm  "  2co  dm 


0) 


The  n’th  row  (#  n  mode)  of  kernel  G,,  is 


rr  ^  -  f  ^,.1! 


Gq  in  matrix  form  is 
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-Pi 


-P2 


£0i  ^  +  SA  ^gpi^ 

^  2coj 

'^02  ^  Qi^  5co2  ^ 


~  rcpi  dn  ^  Qi^  dcp^ 
dm2  2coj  dm2^ 

~  f^02  d/2  ,  Qi^  dco/ 

^n«  I  -V  I 


CO  dmi  2co2  y  ^  5^2  2co„  dmg 


-Pi\? 


r^p^  dVN  1  Qn 

1  ^  (^Off  d/N  1  ^^Pw^ 

CO  dmi  2co^  dmi  j 

1  dm2  2co^  dm2  J 

(#1  mode) 
(#2  mode) 


i#N  mode) 


KNxk) 


Kernel  G„ 


The  last  kernel  is  Ga-  From  (GJ.48),  we  have  the  relationship  between 
Ga  and  Q. 

~  2  Q  ^  ~Q  ^  ]  =  Ga  •  <Jni 

In  a  manner  similar  to  deriving  Gq,  we  can  use  the  Q  and  /  relation  and  the 
partial  derivatives  of  /  to  obtain  the  kernel  Gg  as 

(G  )  =  ^ 

^  2  dm 


gpn  d/  ^  Q-^  dCp^ 
CO  dm  2  dm 


Ga  in  matrix  form  is 


-^1 


-^2 


"gN 


^  Qt  dco,  >1 

~  ffoi  . 

2coj  dmi  J 

a>  dm2 

^  Qt  dco,  ^ 

-  fco,  dr2 

2co2  dmi  j 

fi)  dm2 

1  Qv  dcQj^ 

\  f  (^Off  d/N 

2co^  "dmi  ^ 

j  ®  dm2 

1-1 


'Pi 


'P2 


'On 


'Off 


(#1  mode) 
(#2  mode) 


(#N  mode) 
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Partial  Derivatives  of  c,  Z7,  and  y 


All  foTor  kernels  Gp  Gg  Gq  Gg  are  derived  in  terms  of  partial  derivatives 
of  phase  velocity  c,  group  velocity  U,  and  attenuation  y.  In  this  section,  We 
will  give  all  partial  derivatives  which  will  be  used  in  creating  kernels. 

Perturbation  theory  is  used  to  obtain  phase  velocity  partials  with  respect 
to  medium  parameters,  and  a  numerical  method  which  was  introduced  by 
Rodi  et  al.(1975)  is  used  to  calculate  surface-wave  group-velocity  partial 
derivatives.  The  following  equations  are  calculated  for  a  single  mode  at  a  cer¬ 
tain  frequency,  therefore  we  omit  the  subscript  n  for  n’th  mode.  The  sub¬ 
script  0  is  used  to  denote  the  value  before  introducing  causal  Q.  The  sub¬ 
script  V  can  be  substituted  by  P-wave  velocity  a  or  S-wave  velocity  p.  p  is 
layer  density,  h  is  layer  thickness.  cOr  is  a  reference  angular  frequency  used 
for  introducing  causal  Q. 

Rayleigh  Wave 


1 ,  ^  CO  ^  ^ 

c  =  Cq  +  —  ln( — )  X 


TT  (Or. 


’'=24^ 


^  dco  1 


U  =  Uq\ 


dc  _  3co 
dv  dv 


Cq  Cq 


1 H — —  ln( — ) 

^Qv  COr 
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ap  dp  7C  COr 


dc  dcQ 
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dU 

dh 


dUo  r_^ 
dh  Uo 


Cq  Cq  na 


oh  Co 


on  Co  Cq  ah  nco 


t/p)  ac  2ul  dr 

Co  Co  dQz^  nco  dQ~^ 
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Love  Wave 


1 ,  ,  CO  ,  „  8co 


1-1 


c  =  Co  +  -  ln( — )  X  3^  PQp 
n  (Oj.  op  ^ 


CO  aco 
^  2c2  ^  dyS 


IC  =  M0 


1  +  (2-^^«)(£z£0)+W 

Cq  Cq  JCG) 


8c  8co 


1  1  1  /  ^  \ 
1  +  — pj-ln( — ) 


dc  dco  1  Cl)  ^ 
—  =  3—  +  —  ln( — )| 
op  op  n  (Oj. 


?£0f_A^o-i 

■bp^  2p^^P 


dc  dco 

dc  1  ,  .  CO  .  aco  „ 

^  =  -ln(— )3-y3 

dQjj  n  (Or  dp 

dy  _  acp  1  2y  acp 
dp  2cq  dp  ^  Cp  dp 


CO 


dy 
dp  2cq 


^Cq  (_  P_\Q-i  2/  acp 
8;?  2p  ^  Cp  dp 


dy 

a^ 


Cp  dh 


dy 


CO 


acfi 


dQ/  2cl  dp 
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dU  dUo 


dp  dp 


U  Up  ^c-Cq^  ^2yUo 

C-^o  Cq  Cp  )ro) 


dp  Cp 


2  — -2  — -1 
Cp  C/p 
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—  El  (2-  — ) 

dp  Co  Cq  dp  K(o 


EL 

dp 


ELl  ^  C/q^c-Cq^  ^  2r£/o 

dp  Uq  Cq  Cq  KCO 


+  ^(^^)2 
dp  Co 


E.Ei(2,_ El) ,  ^^0 

dp  Cq  Cq  dp  TtCD 


dU 

dh 


dUQ  U 

dh  Uq 


Uq  (g-cp)  ,  ^yUp 

Cq  Cq  ItCO 


dco 

dh 


(El)2l 

Co 


—  Ei(2-Ei)+ 

dh  Cq  Cq  dh  71(0 


-^0(0  ^0)  ^c  2ul  dy 

Cq  Cq  fTCO  dQ~^ 

4.  Forming  Inversion  Kernel 

As  shown  above,  we  calculate  partial  derivatives  with  respect  to  all 
parameters  (or,  p,  p,  Qp>  h)  instead  of  computing  partial  derivatives  for 
shear  velocity  only.  The  intention  is  that  we  try  to  provide  all  the  possible 
tools  to  interpret  the  seismograms  as  completely  as  possible,  and  as  automat¬ 
ically  as  possible.  Therefore,  it  is  the  users’  responsibility  to  choose  those 
model  parameters  they  want  to  invert;  and  only  the  chosen  part  will  be  used 
to  assemble  the  inversion  kernel. 

To  stabilize  the  inversion,  all  information  for  different  frequencies  inside 
the  inversion  kernel  are  weighted  according  to  their  frequency  amplitudes. 
This  weighting  procedxire  greatly  improves  the  inversion  stability. 

When  inverting  teleseismic  surface-waveforms,  sometimes  the  body 
waves  (e.g.  SS,  SSS)  or  some  unwanted  surface-waves  due  to  improper  rota¬ 
tion  will  have  bad  influence  on  the  cross_correlation  between  the  isolation 
function  and  the  observed  seismogram,  and  cause  signal  misalignment.  To 
avoid  this  problem,  a  window  is  applied  on  the  seismograms  before  doing 
cross-correlation,  and  this  may  successfully  isolate  the  surface  wave 
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wavetrain  from  those  body  waves  outside  the  siuface  wave  wavetrain. 
Although  prewindowing  procedure  is  applied,  the  same  trouble  still  may  hap¬ 
pen  occasionally.  In  such  a  situation,  the  information  at  that  frequency  must 
be  rejected  from  forming  inversion  kernel,  otherwise  it  will  plague  the  inver¬ 
sion  for  its  strong  misaligned  "phase  delay"  or  "group  delay".  Figure  1  shows 
an  example  of  this  problem.  Due  to  improper  rotation,  the  unwanted  Love- 
wave  waveform  appears  on  the  radial  component  seismogram  prior  to  the 
Rayleigh  wave  wavetrain.  This  unwanted  Love  wave  signal  causes  two  kinds 
of  mistakes,  wrong  identification  of  a  Gaussian  wavelet  (Figure  la)  and  sig¬ 
nal  misalignment  (Figure  lb),  which  can  not  be  incorporated  in  inversion. 
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iso 

max=1.55E-03 


Figure  1.  This  is  an  example  showning  how  unexpected  signal  interference  affects  the 
extractmg  procedure  in  GSDF  theory.  There  are  five  traces  presented  to  show  the  different 
processing  stages.  The  top  two  traces  are  the  prefiltered  isolation  filter  and  observed  seismo¬ 
gram,  respectively,  The  tiiird  trace  shows  the  Gaussian  filtered  cross-correlation  at  a  target 
period  of  20.0  seconds.  The  five  extracted  parameters  are  shown.  The  dashed  curves  inside 
the  envelope  are  from  the  S3mthetics  and  the  solid  curve  are  from  the  observed  data.  The 
fourth  trace  is  the  windowed  cross-correlagram  shown  in  the  third  trace.  The  bottom  trace  is 
the  filtered  windowed  cross-correlagram  and  its  five  Gaussian  wavelet  parameters  which  are 
to  be  used  in  further  processing.  Due  to  improper  rotation,  the  Love  wave  appears  on  radial 
component  before  the  Rayleigh  wave  arrival.  The  Love  wave  wavetrain  may  causes  two  Tcind 
troubles  hke  (a)  wrong  determination  of  Gaussian  wavelet  parameters,  and 

We  have  mentioned  that  users  have  the  power  to  decide  what  parame¬ 
ters  will  be  inverted  for  during  the  inversion.  An  aspect  which  strongly 
related  to  this  decision  is  the  phrase  "using  what  kind  information  ?".  During 
each  iteration  in  the  inversion,  we  have  calculated  partial  derivatives  with 
respect  to  parameters  for  each  layer  at  several  appointed  frequencies,  and  all 
this  information  is  rearranged  to  form  four  kind  delays  at  each  frequency. 
After  manually  rejecting  miscalculated  cross_correlations  at  some 
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Figure  1.  (cont’d).  (b)  Misalignment  of  signals  which  produces  significant  bias  phase  and 
group  delay.  This  arises  because  a  period  of  25  seconds  is  considered  compared  to  20  seconds 
in  Figure  la. 

frequencies,  users  have  to  decide  which  combinations  of  6tp,  Stg,  Stq  they 
will  xise  to  invert  for  a  particular  combination  of  model  parameters  (a,  p,  p, 
Qa>  Qfi}  ^)-  From  otm  experience,  we  found  that  when  we  invert  for  velocity 
structure  the  phase  delay  (dtp)  and  the  group  delay  iSfg)  play  the  major  roles 
in  inversion,  and  the  other  two  delays  Stq)  are  better  in  inverting  for 
attenuation  factors.  We  also  found  that  when  the  initial  model  is  far  from  the 
final  result,  using  group  delay  information  in  inversion  can  easily  pull  the 
model  close  enough  to  the  final  model  so  that  phase  delays  can  be  used.  It  is 
only  at  the  final  stage  when  the  S3nithetic  is  very  close  to  the  observed  seis¬ 
mograms  that  the  phase  delay  information  can  be  used  at  a  powerful  fine 
tuning  of  the  model  since  there  is  no  "cycle  skipping"  problem. 
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5.  Synthetic  Test 

A  simple  source  with  strike,  dip,  rake  angles  of  45°,  45°,  45°,  respectively, 
in  20  km  deep  was  used  in  this  synthetic  test.  Receiver  is  located  at  1000  km 
away  along  10°  in  azimuth.  The  "observed"  seismograms  were  created  for  a 
two  layered  crust  with  an  upper  mantle  deviating  slightly  from  the  PREM 
model  (Figure  2).  Here,  we  only  try  to  invert  for  shear  velocity  structure,  so 
all  the  other  parameters  related  to  the  source  and  earth  model  are  assumed 
known.  The  starting  model  is  a  two-layered  crust  model  with  the  PREM 
model  beneath  40  km  depth.  It  is  clearly  seen  in  Figure  3  that  the  starting 
synthetic  seismograms  are  very  far  away  "observed"  seismograms. 


P  P  a  Qb  Qa 


Figure  2.  The  starting  model  (dashed  line)  used  in  synthetic  test  of  GSDF  inversion  algo¬ 
rithm  and  the  "true”  model  (solid  line)  used  to  create  the  "observed  data". 

After  12  iterations,  the  sjmthetic  seismograms  are  almost  the  same  as 
"observed"  seismograms  (Figure  4).  Checking  the  model  differences  between 
the  true  model  and  the  final  model  (Figure  5),  we  can  see  the  crust  struc¬ 
ture  almost  matches  the  "true"  model,  except  in  the  upper  mantle  where  the 
surface  wave  doesn’t  provide  enough  resolving  power  and  the  inverted  struc¬ 
ture  is  causing  the  "zig-zag"  pattern.  The  differences  in  Q  models  are  not  sig¬ 
nificant  at  these  frequencies  and  this  distance. 
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SYNTHETIC  TEST 


START  VEL 


100  150  200  250  300  350  400  450  500 


100  150  200  250  300  350  400  450  500 

FILTER  Flo=0.010  Fhi=0.050  (Hz)  Norder=  4 

LON=0.000  LAT=0.000  DEPTH=  20.00 


AZ=10.000  BAZsIOO.OOO  DIST=  1000.000 


Figure  3.  The  velocity  time  histories  of  "observed  seismograms"  (solid  line)  and  the  those 
from  the  starting  model  (dashed  line)  sure  both  filtered  in  the  frequency  band  0.01-0.05  Hz  by 
using  a  Butterworth  filter  with  four  poles.  The  plotted  seismograms  are  normalized  accord¬ 
ing  the  maximum  amphtude  of  each  component  in  current  frequency  band.  It  is  clear  to  see 
that  the  starting  model  is  not  close  to  the  "true"  model  and  this  may  demonstrate  the  abihty 
of  the  inversion  programs  to  resolve  structure. 

6.  Real  Data  Test 

After  the  successful  synthetic  test,  we  wish  to  test  this  technique  on  real 
data.  The  April  14,  1995  Texas  earthquake  is  the  one  chosen.  According  to 
Sipkin’s  USGS  solution,  the  epicenter  of  the  Texas  earthquake  is  at  30.261°N 
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SYNTHETIC_TEST  FINAL  vel 


LONsO.OOO  LATsO.OOO  DEPTH=  20.00 


AZ=  10.000  BAZ=1 90.000  DISTs  1000.000 


Figure  4.  After  12  iterations,  the  final  inversion  result  show  an  predicted  waveforms  almost 
identical  to  the  "observed  seismograms". 

103.33°W  and  origin  time  is  00:32:55UT.  The  current  collected  data  set  con¬ 
sists  of  49  broadband  stations  from  IRIS,  Canadian  National  Seismological 
Center  (CNSDC),  USGS,  and  UNAM  (Figure  6). 

Some  source  parameters  were  redetermined  on  the  basis  of  fitting  sur¬ 
face  wave  amplitude  spectra.  The  redetermined  source  depth  is  23  km  deep, 
with  strike,  dip,  rake  angles  of  114°,  64°,  -101°,  respectively  with  =  5.6. 
From  some  stations  whose  epicenter  distances  greater  than  3600  km,  it  is 
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p  p  a  Qb  Qa 


Figure  5.  The  comparison  between  the  final  model  and  the  "true  model".  We  can  see  that  the 
2-layer  crust  is  very  close  to  the  "true"  model,  but  that  in  the  upper  mantle  the  model  shows 
some  "zig-zag"  pattern. 

easy  to  see  a  systematically  7-9  seconds  travel  time  difference  between  pP 
and  P  arrivals  (Figure  7).  This  observation  confirms  that  the  source  depth  is 
deeper  than  those  suggested  by  USGS  Sipkin’s  solution  (13  km)  and  Har¬ 
vard’s  CMT  solution  (15  km). 

The  focal  mechanism  generally  matches  the  surface  wave  amplitude 
spectra  radiation  pattern  at  periods  5-60  seconds  (Figures  8a  and  8b).  The 
source  time  function  used  in  this  study  is  a  parabolic  source  time  function 
with  2.0  seconds  duration. 

7.  Inversion  Procedure 

We  have  tried  three  different  inversion  runs,  the  first  two  exhibit  some 
significant  difficulties  in  matching  the  waveform  or  in  reasonableness  of  the 
resulting  model,  so  the  third  one  has  been  adapted  for  inverting  structure.  In 
performing  the  inversions,  we  use  an  earth  flattening  approximation  to  use 
plane-layered  surface-wave  theory  to  generate  synthetics. 

The  first  run  consisted  of  a  joint  inversion  of  both  Rayleigh  and  Love 
wave  seismograms  for  shear  velocity  structure.  The  result  was  that  the 
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210’  220’  230"  240'  250'  260'  270'  280'  290'  300'  310' 

Figure  6.  The  Texas  earthquake  in  relation  to  stations  providing  broadband  data, 
synthetic  seismograms  tended  to  fit  the  largest  surface  wave  amplitude,  and 
ignored  the  small  amplitude  wavefields.  This  arose  because  of  the  amplitude 
level  weighting  used  to  stabilize  the  inversion.  Therefore,  separate  inver¬ 
sions  for  Rayleigh  wave  and  Love  wave  were  necessary. 

The  second  sequence  consisted  of  inverting  for  the  shear  velocity  struc¬ 
ture  from  the  Love  wave,  and  then  using  this  structure  as  an  a  priori  shear 
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Figure  7.  Nine  stations  whose  epicentral  distance  greater  than  3600  km  show  a  general  7-9 
seconds  in  differential  travel  time  between  pP  and  P  phase.  This  supports  a  source  depth 
deeper  than  20  km. 

velocity  structtire  so  the  the  Rayleigh  wave  provided  information  on  the  com- 
pressional  velocity  structure.  This  procedure  can  match  both  Rayleigh  wave 
and  Love  wave  waveforms  very  well,  but  we  found  it  is  impossible  to  find  a 
reasonable  explanation  for  the  anomalous  Poisson’s  ratios  in  the  inverted 
model.  It  is  well  known  that  the  surface  waves  are  insensitive  to  compres- 
sional  wave  velocity  structure,  so  the  only  conclusion  for  this  is  that  the  shear 
velocity  structure  inverted  fi'om  Love  wave  is  not  adequate  for  Rayleigh 
wave.  An  unaccounted  anisotropy  effect  may  cause  the  overcorrection  in  com¬ 
pression  velocity  structure.  Figure  9  compares  waveforms  and  Figure  10, 
shows  the  model  resulting  fi*om  the  second  inversion  procedure.  Although 
the  synthetic  waveform  does  not  perfectly  fit  the  Rayleigh  wave  (Figure  9), 
we  can  see  the  general  featvues  are  matched  using  a  single  model  for  both 
Love  and  Rayleigh  waves.  In  Figure  10  we  see  that  the  P-wave  velocity 
structure  has  lower  values  to  compensate  for  the  high  shear  velocity,  and 
vice  versa.  The  results  are  some  unexplainable  Poisson’s  ratios  such  as  0.135 
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Figure  8.  The  source  mechanism  is  redeterminated  by  fitting  the  surface  wave  amplitude 
spectra  radiation  pattern  at  periods  5-60  seconds,  (a).  The  Love  wave  radiation  pattern  for 
the  preferred  solution.  The  bars  indicate  attenuation  corrected  spectral  amplitudes  in  cm-sec 
normalized  for  geometrical  spreading  to  1000  km. 

for  the  middle  crust  and  0.328  for  the  lower  crust. 

Since  our  forward  S3nithetic  seismogram  algorithm  does  not  include 
anisotropy  effect,  we  modify  our  inversion  procedure  as  as  follows;  we  invert 
^ p_SH  from  Love  wave,  and  another  separate  inversion  of  V p  pgy  from 
Rayleigh  wave  on  the  vertical  component.  We  fixed  the  Poisson^s  ratio  when 
we  invert  V p_psy  from  Rayleigh  wave.  The  fixed  Poisson’s  ratios  are  kept  the 
same  as  the  original  input  model.  In  this  study,  we  set  the  Poisson’s  ratio  as 
0. 25  for  crust,  0. 28  for  the  layers  between  40  km  and  220  km,  and  adopt  the 
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Figure  8.  (Cont’d).  (b)  Comparison  of  observed  and  predicted  Rayleigh  wave  radiation  pat¬ 
terns. 

values  from  PREM  model  for  those  layers  deeper  than  220  km. 

When  we  invert  for  shear-wave  velocity  structures,  the  attenuation  fac¬ 
tors  can  be  determined  simultaneous  either  by  joint  inversion  of  and  Q  or 
by  another  separated  Q  inversion.  The  Q  determination  is  not  definitive, 
since  we  have  a  lot  uncertainties  in  our  source  and  velocity  structures  in  our 
inversion,  so  the  only  objective  criteria  for  determining  Q  is  the  envelope 
shape  of  the  surface  wave  wavetrain.  As  shown  on  Figure  11,  the  Rayleigh 
wave  signal  at  large  distance  is  well  dispersed  with  a  strong  Airy  phase.  This 
Airy  phase  corresponding  is  affected  by  crustal  wave  propagation,  therefore, 
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Figure  9.  The  result  of  inversion  using  the  second  procedure.  This  inversion  procedure 
inverts  S  wave  velocity  from  Love  wave  and  after  that  inverts  P  wave  velocity  from  Rayleigh 
wave  by  assuming  no  amsotropy  effect.  This  shows  a  acceptable  waveform  match  in  phase, 
but  not  in  envelope. 

the  envelope  amplitude  of  the  Airy  phase  is  controlled  by  the  crustal  attenua¬ 
tion  factors.  So  the  Q  structure  determined  is  sensitive  to  the  crust  and 
uppermost  mantle. 
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Figure  10.  The  model  inverted  by  the  second  inversion  procedure. 

8.  Inversion  results 

Data  from  43  of  49  broadband  stations  were  inverted.  In  this  report, 
three  stations  were  selected  to  present  the  inversion  restilts  and  the  wave¬ 
form  fitting  success  will  be  shown  on  several  different  frequency  bands. 

Station  ALE,  which  located  in  Arctic  region  is  the  farthest  station  used 
in  this  study.  There  are  some  interesting  features  for  the  inverted  results. 
Looking  at  the  waveform  fit  in  the  low  frequency  range  (0.005-0.03  Hz  band¬ 
pass;  shown  as  Figure  11a),  it  is  clear  to  see  that  the  S3mthetic  seismograms 
successfully  match  the  observed  surface  wave.  Two  velocity  models,  PSV  and 
SH  are  inverted  respectively  from  Rayleigh  wave  and  Love  wave. 

For  ALE,  the  attenuation  factor  was  fixed  dining  inversion.  A  high  Q 
structure  (Qs  =  1000)  is  used  for  crust,  a  low  Q  (Qg  =  100)  was  adopted  for 
structvire  between  40  km  and  500  km,  and  using  Qg  =  143  for  those  deeper 
than  220  km.  And  from  Figures  lla,b,  the  synthetic  Rayleigh  wave  and  Love 
wave  amplitudes  only  show  small  deviation  from  observed  seismogram, 
therefore  the  Q  model  is  considered  adequate.  We  also  notice  this  sharp  Q 
contrast  between  crust  and  mantle  is  a  common  character  for  those  stations 
located  inside  the  North  America  craton  region. 
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Figure  11a.  Waveform  fit  of  the  final  inverted  model  for  ALE  in  the  frequency  band  of 
0.005-0.03  Hz.  The  signal  which  arrives  at  1020  seconds  and  around  1300  seconds  are  S  and 
SS  phase,  respectively.  From  the  S  phase  waveform,  we  can  say  that  the  source  time  func¬ 
tion  used  in  this  inversion  is  a  little  short  but  is  close  enough. 

The  second  station  is  FCC  which  located  on  the  west  shoreline  of  Hudson 
Bay  and  is  in  the  center  of  the  North  America  craton.  The  inverted  models 
(Figure  13)  show  high  shear  velocity  for  upper  mantle  but  not  as  high  as  SNA 
model  (Grand  and  Helmberger,  1984).  The  current  inversion  result  shows 
4.5%  anisotropy  between  70  km  and  140  km.  The  synthetics  fit  data  well  at 
low  frequency  (Figure  14a)  and  can  fit  the  fundamental  mode  as  high  as  0.1 

Hz  (Figure  14  ),  but  have  difficulty  to  produce  some  higher  mode  arrivals  at 
800  and  940.  the  time  range. 

The  final  station  is  PAS.  The  wave  propagates  through  the  southern 
Rocky  Mountains.  The  inverted  model  (Figure  15)  does  not  require 
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Figure  11b.  (Cont’d).  Waveform  fit  in  the  0.005-0.05  Hz  frequency  band, 
anisotropy  and  the  shear  velocity  model  is  very  close  to  the  TNA  model 
(Grand  and  Helmberger,  1984).  The  synthetics  fit  the  S  phase  which  arrives 
around  350  seconds  and  surface  wave  very  well  (Figures  16abc)  and  this  sug¬ 
gests  that  there  is  a  velocity  discontinuity  at  220  km  which  can  not  be  seen  in 
the  TNA  model.  The  Q  is  very  low,  with  the  average  Q  for  crust  of  lower  than 
200.  However  we  found  one  interesting  feature  about  the  Q  behavior.  The  Q 
values  between  40  and  220  km  are  slightly  higher  than  PREM  model,  i.e.,  a 
low  Q  crust  underlain  a  slightly  high  Q  upper  mantle  with  respect  to  the  ref- 
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Figure  12.  The  inverted  models  for  ALE.  TWo  models  (PSV  and  SH)  were  inverted  for  Love 

there  is  one  slight  anisotropic  zone  above 


P  P  ot  Qb  Qa 


“verted  model  for  FCC.  There  is  4.5%  anisotropy  effect  exists  between 

/U  and  140  inn. 


erence  model. 


TEXAS950414 


F  C  C  VEL 


600  700  800  900  1000  1100  1200 

I  '  '  '  '  I  ■  '  '  '  I  '  '  ■  '  I  '  '  '  '  I . *-*-1 


LON—1 03.327  LAT-30.261  DEPTH-  23.00 

AZ-9.879  BAZ-1 96.549  DIST-  3244.330 

Figure  14.  The  waveform  fitting  for  final  model  in  FCC  are  shown  at  three  frequency  bands  : 
(a)  0.01-0.03  (b)  0.01-0.05  (c)  0.01-0.1  Hz.  The  SS  phase  arrives  at  660  seconds.  The  inverted 
model  can  fit  fundamental  mode  Love  wave  and  Rayleigh  wave  waveforms  as  high  as  0.1  Hz, 
but  it  lacks  the  ability  to  simulate  the  higher  modes. 

9.  Conclusion 

We  implemented  the  Generalized  Seismological  Data  Functionals  tech¬ 
nique  of  Gee  and  Jordan  (1992)  in  a  surface-wave  waveform  inversion  algo¬ 
rithm.  A  simple  S3nathetic  test  shows  its  robust  inversion  ability.  After  a  suc¬ 
cessful  S3mthetic  test,  we  useed  this  inversion  algorithm  on  real  data  to  fur¬ 
ther  test  its  ability.  The  Texas  earthquake  (30.26  °N  103.33°W,  00:32:55UT, 
April  14,  1995)  is  a  very  good  earthquake  for  this  purpose  because  it  was  well 
recorded,  its  source  depth  is  well  constrained,  the  focal  mechanism  is  fairly 
determined,  and  the  seismic  moment  is  constrained  by  long  period  surface 
waves.  The  inversion  results  are  excellent  and  show  some  interesting  fea¬ 
tures.  For  the  craton  there  is  some  evidence  for  anisotropy  and  crustal  Q  is 
high.  For  the  mountain  region,  although  the  inverted  model  show  very  simi¬ 
lar  shear  velocity  structure  like  TNA  model,  but  from  the  S-wave  waveform 
the  model  prefers  a  velocity  discontinuity  at  220  km.  These  features  worth 
more  effect  on  them. 
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Figure  14.  (Cont’d).  (b)  0.01-0.05  Hz. 
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Figure  14.  (Confd).  (c)  0.01-0.1  Hz.  The  SS  phase  arrives  at  660  seconds.  The  inverted 
model  can  fit  fundamental  mode  Love  wave  and  Rayleigh  wave  waveforms  as  high  as  0.1  Hz, 
but  it  lacks  the  ability  to  simulate  the  higher  modes. 
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Figure  15.  The  inverted  models  for  PAS.  There  is  no  clear  anisotropy  effects.  The  crustal  Q 
is  very  low. 
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Figure  16.  Waveform  fitting  are  shown  on  three  frequency  bands  :  (a)  0.01-0.03  (b)  0.01-0.05 
and  (c)  0.01-0.1  Hz.  The  S  wave  signal  arrives  at  350  seconds. 
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Figure  16.  (Cont’d).  (b)  0.01-0.05  Hz. 
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Figure  16.  (Cont’d).  (c)  0.01-0.1  Hz. 
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